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ABSTRACT 

This  thesis  analyzes  the  sway,  yaw,  and  roll  dynamic  stability  of  neutrally 
buoyant  submersibles.  Utilizing  the  hydrodynamic  coefficients  for  a  Mark 
IX  Swimmer  Delivery  Vehicle  (SDV)  as  a  base-line  model,  the  linearized 
equations  of  motion  for  the  decoupled  steering  and  roll  systems  are 
compared  to  the  coupled  system.  Two  different  configurations  of 
hydrodynamic  coefficients  are  considered  along  with  the  effects  of  varying 
the  vertical  (Zg)  and  longitudinal  (Xg)  centers  of  gravity  of  the  vehicle  while 
the  longitudinal  center  of  buoyancy  (Xb)  is  held  constant.  Results 
demonstrate  the  significant  effects  on  stability  of  coupling  the  steering  and 
roll  equations  of  motion,  and  the  importance  of  Zg  and  Xg  selection  in 
minimizing  those  effects  while  retaining  stability.  Perturbation  analysis 
results  confirm  the  essential  dependence  of  the  linearized  coupled 
equations  on  the  separation  of  Xg  and  Xb- 
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I.    INTRODUCTION 
A.     GENERAL 

As  the  missions  for  submersibles  become  more  complex  and 
demanding,  the  requirement  for  a  highly  stable  platform  becomes 
increasingly  important  so  the  operator(s)  can  concentrate  on  matters  other 
than  station-keeping.  Submersible  simulators  have  not  been  employed  to 
any  great  extent  thus  far,  since  the  actual  system  is  relatively  inexpensive 
and  the  dynamics  are  usually  very  non-linear  and  difficult  to  model.  The 
analysis  of  a  submersible  can  be  significantly  more  complex  than  the 
analysis  of  a  conventional  submarine  or  aircraft,  since  the  presence  of 
ancillary  equipment  such  as  manipulators,  video  devices,  and  tethers 
introduce  extra  cross-coupling  terms  usually  absent  in  other,  more 
symmetric,  vehicle  dynamic  analyses.  [Ref.  1]  Additionally,  all 
mathematical  models  include  simplifying  assumptions  and  errors  in  the 
models  hydrodynamic  coefficients. 

Submersibles  typically  have  a  variety  of  complex  dynamic  interactions 
that  can  severely  inhibit  maneuverability  and  control  performance.  The 
goal  of  this  thesis  is  to  present  an  understanding  of  the  coupling  effects  on 
straight  lira  motion  stability  in  the  horizontal  plane  using  a  linearized 
model,  and  the  primary  means  of  minimizing  these  effects.  Development  of 
the  mathematical  models  for  both  the  coupled  and  uncoupled  maneuvering 
and  roll  equations  of  motion  is  presented  in  Chapter  IT.  Utilizing  two 
different  configurations  of  hydrodynamic  coirfTiripmt.s,  the  degree  of  stability, 
regions  of  stability,  and  linear  simulations  for  the  coupled  and  uncoupled 


systems  are  presented  in  Chapter  III.  A  Hamming  method  nonlinear 
simulation  [Ref.  2],  which  is  similar  to  a  fourth  order  Runge-Kutta 
integration  technique,  is  conducted  on  both  configurations  to  compare  with 
the  results  obtained  for  the  linearized  models.  Chapter  IV  develops  a 
perturbation  analysis  to  demonstrate  the  strong  degree  of  dependence  on 
the  separation  between  the  longitudinal  centers  of  buoyancy  and  gravity  to 
the  solution  of  the  linearized  coupled  equations  of  motion.  Chapter  V 
summarizes  the  results  and  provides  recommendations  for  future 
submersible  modelling  research.  Appendix  A  contains  the  computer 
programs  utilized  for  the  linear  and  nonlinear  simulations. 

B.    PARAMETER  DEFINITION 

The  values  for  the  hydrodynamic  derivatives  and  vehicle  dimensions 
are  from  Smith,  Crane,  and  Summey  [Ref.  3],  with  the  following  exceptions: 

•  Yr  -  the  force  in  sway  due  to  a  change  in  yaw  rate 

•  Nv  -  the  moment  due  to  a  change  in  sway  velocity. 

These  two  coefficients  were  modified  to  produce  two  different  models  that 
would  have  one  eigenvalue  change  sign  for  a  reasonable  range  of 
longitudinal  and  vertical  centers  of  gravity.  A  comparison  between  the 
actual  non-dimensional  values  and  those  used  in  this  thesis  is  as  follows: 


X. 


Configuration  'A' 
Configuration  *B' 
Actual  SDV 


-3.500E-02 
-5.940E-02 
2.970E-02 


Nv_ 

-1.484E-03 
-1.484E-03 
-7.420E-03 


The  effects  of  changing  these  coefficients  are  illustrated  in  Figures  1  and  2 
on  the  following  pages.  Additionally,  the  analysis  presented  herein  is 
conducted  in  dimensional  form;  hence  the  nominal  forward  longitudinal 
velocity  *U'  appears  in  the  equations  of  motion.  All  calculations  and 
simulations  utilize  a  value  of  5  ft/sec  for  'U'.  The  coordinate  system 
convention  is  the  standard  body-fixed,  right  -hand  orthogonal  axis  system 
employing  the  Euler  angle  approach. 

1.   Variables 

x,  y,  z  Distances  along  the  body  fixed  principal  axes. 

u,  v,  w  Translational  velocity  components  of  model  relative  to 

fluid  along  body  axes. 

p,  q,  r  Rotational  velocity  components  of  model  along  body  axes. 

X,  Y,  Z  Hydrodynamic  force  components  along  body  axes. 

K,  M,  N  Hydrodynamic  moment  components  along  body  axes. 

\j/,  6,  (j)  Yaw,  pitch,  and  roll  angles;  positive  values  following  the 

right-hand  rule. 

Xg,  Yg,  Zg  Center  of  gravity  coordinates  along  principal  axes. 

Xb,  Yb,  Zb  Center  of  buoyancy  coordinates  along  principal  axes. 

Ixx.  Iyy>  Izz         Moments  of  inertia  about  the  principal  axes. 

XnoSe»  Xtaii  Distances  from  body  center  along  the  longitudinal  axis 

used  in  the  crossflow  force  and  moment  integrals. 
Values   are   located   within   the   nonlinear   computer 
simulation  program  in  Appendix  A. 

h(x),  b(x)  Model  width  and  height  values  used  in  the  crossflow 

force  and  moment  integrals.   Values  are  located 
within  the  nonlinear  computer  simulation  program 
in  Appendix  A. 
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II.   STABILITY  OF  MOTION 

A.    EQUATIONS  OF  MOTION 

The  horizontal  plane,  nonlinear  equations  of  motion  for  a  submersible 
as  developed  by  Smith,  Crane,  and  Summey  [Ref.  1]  are  shown  below  in 
Equations  (2.1). 

Sway:  m[v  +  ur-  wp  +  Xg(pq  +  r)-  Yg(p2  +  r2)  +  Zg(qr-p)]  = 

(2.1a) 

[Ypp  +  Yrr  +  Ypqpq  +  Yqrqr]  +[Yvuv  + Yvwvw  + Y5ru25r]    + 

[Yvv  +  Ypup+  Yrur+  Yvqvq  +  Ywpwp  +  Ywrwr]  +    (W  -  B)cos  6  sin  0  - 

J"  [CDyh(x)(v  +  xr)2  +  CD2b(x)(w-xq)2l(v  +  xr)dx 

Ucf(x) 


Yaw:  lat  +  (Iyy-Ixx)pq  -  I^p2^2)  -  IyZ(pr+q)  + 


(2.1b) 


Ixz(qr-p)  +  m[Xg(v  +  ur-wp)  -  Yg(u-vr+wq)]  = 
[Npp  +  Nrf  +  Npqpq  +  Nqrqr]  + 
[NyV  +  Npup  +  Nrur  +  Nvqvq  +  Nwpwp  +Nwrwrj  + 
[Nvuv  +  Nvwvw  +  Ngj-u^rj  +  (XgW-XbB)cos  0  sind)  +  (YgW-YbB)sinO 


J  [CDyh(x)(v  +  xr)2   +  CDzb(x)(w-xq)2](v+xr)  (x)  dx  +  u2N 


Ucf<x) 


prop 


Roll:  IxxP+qrdzz  -Iyy)  +  IXy(pr-q)-Iyz(q   -r2)- 

(2.1c) 

Ixz(pq  +r)  +  m[Yg(w-uq  + vp)-Zg(v  +  ur-wp)]   = 

[Kpp  +  Krf  +  Kpqpq  +  Kqrqr]  + 

[Kvv  +Kpup  +  Krur+  Kvqvq  +  Kwpwp+Kwrwr]  + 

[Kvuv+Kvwvw]  +  (YgW-YbB)cos6cos(t)  - 

(ZgW-ZbB)cos  6  sin  <j)  +u2Kprop 


B.    SIMPLIFICATIONS 

In  order  to  obtain  the  linearized  equations  of  motion  about  a  level  flight 
path,  the  following  simplifications  were  utilized: 


•  The  translational  velocity  (w)  and  acceleration  (w  )  in  the  z-direction 
are  zero. 

•  The  rotational  velocity  (q)  and  acceleration  (q)  in  the  y-direction  are 
zero. 

•  The  acceleration  in  the  longitudinal  direction  (li  )  is  zero. 

•  The  cross-products  of  inertia  are  zero  by  virtue  of  a  body-centered 
coordinate  system. 

•  The  submersible  is  neutrally  buoyant  so  B  =  W. 

•  The  longitudinal  center  of  buoyancy  (Xb)  and  the  vertical  center  of 
buoyancy  (Zb)  are  located  at  the  origin  of  the  body-fixed  coordinate 
system. 

•  The  lateral  center  of  buoyancy  (Yb)  and  the  lateral  center  of  gravity 
(Yg)  are  located  at  the  origin  of  the  body-fixed  coordinate  system. 

•  Dynamic  stability  analysis  is  performed  with  all  controls  fixed; 
hence,  all  forces  and  moments  due  to  control  surfaces  are  zero. 

•  The  angle  of  pitch  (6)  is  sufficiently  small  for  sin(0)  to  equal  zero. 


•     From  Smith,  Crane,  and  Summey  [Ref.  1]  the  propeller  coefficients 
Kprop  and  Npr0p  are  zero. 


C.   COUPLED  STABILITY  EQUATIONS 

When  the  simplifying  assumptions  from  the  preceeding  section  are 
applied,  the  resulting  linearized  equations  are: 

(2.2a)  Sway:  m[v +ur  +  Xg(r)-Zg(p)]    =  Yf 

(2.2b)  Yaw:  I^r  +  mXg(v+ur)   =  Nf 

(2.2c)  Roll:  Ixxp  -    mZg(v+ur)    =  Kf 

where  the  force  and  moment  representations  are  given  by: 
(2.3a)  Sway  Force:        Yf   =  Yvv  +Yvv  +  Yj.r  +  Yrr  +  Ypp  +  Ypp 

(2.3b)  Yaw  Moment:  Nf   =  Nvv +Nvv  +  Nrf +  Nrr  + Npp  +  Npp  + (XgW-XbB)(p 

(2.3c)  Roll  Moment:     Kf   =  K^v  +  Kvv  +  K^r  +  Krr  +  Kpp  +  Kpp+  (ZgW-ZbB)cp 

Equations    2.2    and    2.3    may   be    combined   to   form   the    state    space 
representation: 

Ax  =  Bx  (2.4a) 

where  x  =  [p    0   v     r]  (2.4b) 


Ixx-Kp  0    -(Kv  +  MZg)  -Kr 

0  10  0 

-(Yp+MZg)  0        M-Yv  MXg-Yr 

-NP  0      MXg-Nv  Izz-NrJ 


(2.4c) 


KPU     ZbB-ZgW  KvU  U(MZg+Kr) 
10               0  0 

YpU  0  YvU  U(Yr-M) 

NPU    XgW-XbB  NvU  U(Nr-MXg). 


(2.4d) 


The  stability  of  the  coupled,  linearized  system  depends  on  the  location  of  the 
four  eigenvalues  of  det(B  -  XA)  =  0,  which  has  a  characteristic  equation  of 
the  form 

A?i4  +  BX3  +  C?i2  +  D?i  +  E  =  0,  (2.4e) 

where  the  coefficients  are  complicated  permutations  of  the  elements  in 
matrices  A  and  B.  The  values  for  A,  B,  C,  D,  and  E  are  given  in  Equations 
(2.4f  -  2.4j)  in  terms  of  lower  case  letters  that  represent  entries  in  matrices  A 
and  B;  they  are  explicitly  defined  in  Table  1. 


A  =  a(j-u  -  1-r)  +  d-(u-h  +  o-l)  -  f-(r-h  +  o-j)  (2.4f) 

B  =  e-(uh  +  o-l)  -  d-(u-i  -  w-h  +  o-m  -  p-1)  -  a-(j-w  +  k-u  -  1-x  -  r-m) 

-  b-(j-u  -  1-r)  +  f-(r-i  -  x-h  +  o-k  -  p-j)  -  g-(r-h  +  o-j)  (2.4g) 

C  =  a(kw  -  mx)  +  b(jw  +  k-u  -  1-x  -  r-m)  +  g(ri  -  x-h  +  ok  -  pj) 
+  c(ju  -  1-r)  -  d-(q-l  -  w-i  +  m-p)  -  e(ui  -  w-h  +  om  -  p-1) 
+  f-(q-j  -  x-i  +  p-k)  (2.4h) 


D  =  g-(qj  -  x-i  +  pk)  -  f-q-k  +  d-q-m  -  e-(ql  -  wi  +  mp)  -  b(k-w  -  mx) 

-  c-(j-w  +  ku  -  lx  -  r-m)  (2.4i) 

E=    c(kw  -  mx)  +  e-q-m  -  g-q-k  (2.4j) 


TABLE  1 

COEFFICIENT 

DESCRIPTOR  VALUES 

a  =   Lcx-Kp 

b  =  KPU 

c  =  ZgW-ZbB 

d  =  MZg+Kv 

e  =  KvU 

f  =  Kf 

g  =  MZgU 

h  =  MZg  +  Yp 

i  =  YPU 

j   =  M-Yv 

k  =  YvU 

1    =    MXg-Yr 

m  =  U(Yr-M) 

o  =  NP 

p  =  NPU 

q  =  XbB-XgW 

r  =   MXg-Nv 

U    =    Izz  -  Nr 

W    =    U(Nr-MXg) 

x  =  NvU 

D.   UNCOUPLED  STABILITY  EQUATIONS 

Using  matrices  A  and  B  from  the  preceeding  section  (Equations  2.4c  and 
2.4d),  uncoupling  the  steering  and  roll  equations  is  straightforward: 


STEERING  EQUATIONS 


M-Yv       MXg-Y/ 
MXg-Nv      Ia-Nt. 


YvU       U(Yr-M) 
LNvU    U(Nr-MXg)JLrJ 


(2.5a) 


ROLL  EQUATIONS 


Ixx-Kp    0 
0         1 


KpU    ZbB-ZgW 

1  0 


P 

L<!>J 


(2.5b) 
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The  characteristic  equations  for  the  uncoupled  steering  and  roll  equations 
are  much  simpler  than  that  for  the  coupled  system,  and  are  given  as 
Equations  (2.5c  and  2.5d)  in  terms  of  the  hydrodynamic  parameters  rather 
than  the  descriptors  used  in  the  previous  section. 

The  steering  characteristic  equation  form  is:  Al^2  +  BlA.  +  Cl  =  0,  where 
AL   =  (M-Yv)(Izz-Nr)  -  (MXg-Yr)(MXg-Nv) 

BL   =   -[azz-Nr)(YvU)+(M-Yv)(Nr-MXg)(U)+  (2.5c) 

(Yr  -MXNv  -MXg)(U)+(Yr  -MXg)(NvU)] 


CL  =  (YvU2)(Nr-MXg)  -  (NvU2)(Yr-M) 
The  roll  characteristic  equation  form  is:   ArA.    +  Br?i  +  Cr  =  0,  where 


Ar 

=  Lex  -Kp 

Br 

=    -KpU 

Cr 

=    zgw 

(2.5d) 
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III.  RESULTS  AND  DISCUSSION 

A.  DEGREE  OF  STABILITY 

The  effects  of  changing  the  longitudinal  and  vertical  centers  of  gravity 
(while  keeping  the  center  of  buoyancy  fixed  at  the  vehicle  center)  on  the 
degree  of  stability  for  configuration  'A'  is  illustrated  in  Figure  3.  Degree  of 
stability  as  utilized  in  this  thesis  is  defined  as  the  maximum  real  value  of  all 
characteristic  roots,  with  negative  values  indicating  a  stable  situation.  The 
degree  of  stability  for  the  uncoupled  system  is  represented  by  the  dashed  line. 
It  can  be  seen  that  the  critical  value  of  Xg  for  which  motions  become  unstable 
is  clearly  a  function  of  the  metacentric  height  Zg,  whereas  the  uncoupled 
system  predicts  a  constant  value  of  Xg.  Figure  4  displays  the  variation  of  the 
imaginary  part  of  the  root  value  for  configuration  'A'. 

Figures  5  and  6  are  analogous  to  Figures  3  and  4  for  configuration  'B'; 
they  show  the  degree  of  stability  for  varying  Xg  and  Zg  values.  For  this 
configuration,  the  degree  of  stability  has  a  stronger  dependence  on  the 
location  of  Xg  and  the  value  of  Zg.  For  almost  all  positive  values  of  Xg,  the 
complex  conjugate  roots  are  increasing  in  value  and  eventually  becoming 
positive;  this  indicates  an  oscillatory  response  that  diverges  when  the  degree 
of  stability  is  positive. 

B.  REGIONS  OF  STABILITY 

Figure  7  shows  the  region  of  stability  for  configuration  'A',  with  the 
uncoupled  system  represented  by  a  dashed  line.  The  uncoupled  system 
predicts  stability  for  all  values  of  Xg  greater  than  0.18  ft,  while  the  coupled 
system  predicts  an  additional  region  enclosed  by  the  triangular  area  to  the 
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left  of  the  dashed  line.  Figure  8  is  analogous  to  Figure  7  for  configuration  'B'. 
The  large  discrepancies  between  the  predicted  regions  of  stability  in  this  case 
occur  for  Zg  values  less  than  0.045  ft.  For  small  values  of  Zg,  there  are 
corresponding  small  regions  of  Xg  where  stability  is  predicted  by  the  coupled 
equations  but  not  by  the  uncoupled  equations.  The  root  values  in  this  region 
are  complex  conjugates  with  very  small  real  parts.  Figures  9(a)  and  9(b) 
illustrate  the  effect  of  co-locating  Xb  and  Xg.  This  results  in  a  degree  of 
stability  for  the  coupled  system  that  is  nearly  identical  to  the  uncoupled 
system. 
C.  INTERPRETATION  OF  RESULTS 

It  may  be  shown  by  applying  Routh's  criterion  [Ref.  4:  pp.  211-218]  that 
for  a  fourth  order  equation  of  the  form  A?i4  +  BX  +  C^.2  +  DX  +  E  to  have 
roots  with  all  negative  real  parts  the  following  must  apply: 

i.)  BCD  -  AD2  -  EB2  >  0 
ii.)  E  >  0. 
If  the  quantity  'E'  is  less  than  zero,  the  system  will  become  unstable  and  the 
resulting  motion  will  be  a  simple  divergence.  If,  however,  the  value  of  the 
quantity  BCD  -  AD2  -  EB2  is  negative,  the  resulting  instability  will  result  in 
an  oscillatory  motion  due  to  the  presence  of  complex  conjugate  roots  with 
positive  real  parts. 

For  the  coupled  system  of  equations,  the  condition  E  =  0  yields: 
Zg  =  Xg  .[Kv(YrM)/(YvNr  -  Nv(Yr-M))]. 
while  the  uncoupled  system  of  equations  reduces  to  a  constant  term 
expression  for  Xg: 

Xg  =  [(MYv)/(YvNr  -  Nv(YrM))]. 
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This  explains  the  differences  in  the  regions  of  stability  illustrated  in  Figure  7, 
since  the  above  equations  show  a  linear  relationship  between  Zg  and  Xg  for 
the  coupled  equations  and  a  constant  value  for  Xg  for  the  uncoupled 
equations.  When  the  value  for  Xg  coincides  with  the  value  for  Xt>,  the 
constant  term  'E'  in  the  coupled  equations  is  reduced  to  that  of  the  constant 
associated  with  the  uncoupled  equations,  and  the  resulting  predicted  degree 
of  stability  no  longer  depends  on  the  value  of  Zg.  Substituting  the  coefficient 
values  for  'E'  serves  to  clearly  illustrate  the  reduction: 

E    =      (ZgW)[(YvU2)(NrMXg)  -  (NvU2)(Yr-M)]     + 
(KvU2)(Yr-M)(XbB-XgW)  -  (MZgUXYvUXXbB-XgW). 

When  Xb  =  Xg  and  B  =  W  for  the  neutrally  buoyant  case,  the  second  and  third 
terms  are  reduced  to  zero.  When  'E'  is  then  set  equal  to  zero  (the  condition 
for  determining  where  the  real  roots  change  from  negative  to  positive),  the 
dependence  on  the  value  Zg  is  removed  and  the  expression  for  Xg  is  the 
constant  given  for  the  uncoupled  equations.  This  reduction  is  the  explanation 
for  the  appearance  of  Figures  9(a)  and  9(b).  When  the  longitudinal  centers  of 
buoyancy  and  gravity  coincide  for  a  neutrally  buoyant  vehicle,  the  degree  of 
stability  for  the  coupled  system  of  equations  covering  all  metacentric  height 
values  is  equivalent  to  the  uncoupled  system  of  equations. 

A  simple  reduction  of  the  equation  resulting  from  Routh's  criterion  to 
determine  when  a  pair  of  complex  conjugate  roots  crosses  the  zero  axis  is  not 
easily  accomplished.  Figure  10  is  presented  as  confirmation  that  the 
locations  of  Xg  for  which  BCD  -  AD2  -  EB2  =  0  matches  the  locations  given 
graphically  in  Figure  5. 

22 


cc 


5 

c 

o 
U 


be 


tc 
X 

CO 

> 

c 


tevovMazvHMcDH)  aaznvraoN 


23 


D.  LINEAR  SIMULATIONS 

Figures  11  through  14  present  comparisons  of  the  coupled  and  uncoupled 
system  responses  for  configuration  'A'.  For  the  stable  cases  Xg  =  0.40  ft  and 
for  the  unstable  cases  Xg  =  -0.20  ft,  while  Zg  is  0.20  ft  for  both.  The  unstable 
coupled  case  (Figure  11)  illustrates  a  simple  divergence  for  both  angle  of  drift 
and  roll  angle.  This  is  expected  since  the  roots  are  not  complex  conjugates. 
The  unstable,  uncoupled  case  accurately  predicts  the  divergence  in  angle  of 
drift,  but  the  roll  response  is  predicted  to  be  stable.  This  may  be  explained  by 
examining  the  uncoupled  equation  of  motion  in  roll: 

(J)'*  -  0'  (KpU)/(Ixx  -Kp)  +    0  (ZgW)/(Ixx  -Kp)  =  o. 
Substituting  values  for  configuration  'A'  yields: 

0"  +  (J)' (1.475)  +  (J)  (0.720)  =  0. 
The  solution  of  this  equation  results  in  a  natural  frequency  of  0.85  rad/sec 
and  a  damping  factor  of  0.869.  This  is  an  underdamped  case,  since  both  roots 
have  negative  real  parts  and  are  complex  conjugates.  Substituting  values  for 
configuration  'B'  yields: 

0"  +  0' (1.475)  +  0  (0.144)  =  0. 
The  solution  for  this  case  results  in  a  natural  frequency  of  0.380  rad/sec  and  a 
damping  factor  of  1.941,  which  represents  an  overdamped  situation.  This 
also  demonstrates  that  a  vehicle's  natural  frequency  in  roll  may  be  increased 
by  increasing  the  metacentric  height.  The  effects  of  the  underdamping  may 
be  seen  in  Figures  12  through  14. 

Figures  15  through  20  provide  a  comparison  between  the  coupled  and 
uncoupled  equations  for  configuration  'B'.  The  effects  of  the  overdamping  are 
evident  in  figures  15  through  17.     This  set  of  figures  demonstrate  the 
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Xg 

0.20 

0.20 

Roll  Stable 

Y 

Y 

Drift  Stable 

Y 

N 

magnitude  of  the  discrepancies  when  the  coupling  effects  are  not  considered. 
The  comparison  below  summarizes  the  results  of  the  figures  representing 
simulations  for  configuration  'B\  where  'C  stands  for  coupled  and  'UC'  for 
uncoupled. 

Q UC Q UC Q H£ 

1.00  1.00  1.50  1.50 

N  Y  N  Y 

N  N  N  Y 

As  was  seen  before,  the  effects  of  the  coupling  on  the  system  results  in 
complex  conjugate  eigenvalues  with  positive  real  parts.  Therefore,  the  larger 
values  of  Xg  result  in  increasingly  divergent  oscillations  instead  of  the 
stability  predicted  by  the  uncoupled  system. 

Figure  21  is  a  three-dimensional  presentation  of  the  roll  amplitude  vs 
time  for  configuration  'B'  as  Xg  varies  from  0.15  ft  to  1.50  ft.  This  mesh 
capability  of  MATLAB  allows  a  comparative  view  of  several  solutions,  and 
the  behavior  of  the  roll  response  for  the  coupled  system  is  easier  to  discern. 

E.  NON-LINEAR  SIMULATIONS 

In  order  to  provide  a  measure  of  the  accuracy  of  the  results  obtained 
utilizing  the  linearized  equations  of  motion,  a  simulation  program  for  the 
non-linear  equations  of  motion  was  developed  using  Hamming's  method 
[Ref.  2].  Hamming's  method  utilizes  a  Milne  predictor  and  incorporates  a 
modifier  step  prior  to  the  correction  step.  The  primary  advantage  of  using 
Hamming's  method  is  that  only  two  derivative  function  evaluations  are 
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required  per  step  rather  than  the  four  or  more  normally  required  by  other 
popular  methods.  The  local  error  is  of  the  same  order  of  magnitude  (h4)  as  a 
more  time  consuming  process  such  as  a  fourth  order  Runge-Kutta,  but  the 
reduction  in  function  evaluations  results  in  a  faster  simulation.  The  formula 
is  presented  below,  and  may  also  be  found  in  the  nonlinear  computer 
program  simulation  in  Appendix  A. 

HAMMING'S  METHOD 

yG+Dpredicted  =  y(i-3)  +  (4h/3)[2fli)  -  fli-1)  +  2fli-2)] 

y(i+l)modified  =   y(i+Dpredicted  -  (1 12/12 l)[y(i)predicted  -  jKi Corrected] 
ytf+Dcorrected    =    d/8)[9y(i)  -  y(i-2)  +  3h{fU+l Modified  +  2f(i)  -  fli-l)}] 
y(i+l)    =   Vft+Dcorrected   +   (9/121)[y(i+l)predjcted  "   y(i+Dcorrected] 

The  first  four  values  must  be  determined  by  another  method;  the  Euler 
linear  solution  with  a  small  step  size  *h'  proved  sufficient. 

Figure  22  represents  the  simulation  for  the  unstable  representation  of 
configuration  'A'.  Rather  than  the  exponentially  increasing  roll  angle  and  the 
-90  angle  of  drift  computed  with  the  linear  simulation,  the  nonlinear  solution 
predicts  an  angle  of  drift  that  reaches  -15  degrees,  and  then  slowly  diverges. 
The  roll  angle  reaches  a  steady  state  value  of  approximately  three  degrees. 
Figures  23  and  24  are  the  nonlinear  simulation  results  for  stable 
configurations  of  'A'  and  'B',  respectively.  They  are  nearly  identical  to  the 
results  obtained  using  the  linear  simulation  and  displayed  as  Figures  12  and 
18. 

Figure  25  is  the  simulation  for  an  unstable  configuration  *B\  Quite 
notable  are  the  steady  state  roll  angle  and  angle  of  drift  after  approximately 
250  seconds  rather  than  the  exponentially  increasing  divergence  apparent  in 
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the  corresponding  linear  simulation.  Figure  26  provides  a  comparison 
between  the  linear  and  nonlinear  solution  for  angle  of  roll.  The  enlarged 
view  of  the  steady  state  region  provided  in  Figure  27  better  illustrates  the 
constant  limits  of  oscillation.  By  plotting  the  roll  angle  versus  the  angle  of 
drift,  an  elliptical  trajectory  is  apparent  in  Figure  28.  This  result  is  similar 
to  the  results  obtained  by  Schmidt  and  Wright  in  their  analysis  of  high 
performance  aircraft  wing-rock  [Ref.  5].  They  postulate  that  a  possible 
explanation  for  the  limit  cycle  is  the  inertial  coupling  between  a  stable 
longitudinal  and  an  unstable  lateral  mode.  Similar  results  in  the  work  are 
attributed  to  dynamic  as  well  as  hydrodynamic  coupling.  Figures  25  and  28 
indicate  that  the  nonlinear  interactions  for  the  unstable  conditions  of 
configuration  'B'  provide  a  significant  amount  of  damping  to  the  rolling 
motion.  The  limiting  of  the  rolling  motion  accomplished  by  including  the 
nonlinear  terms  then  serves  to  limit  the  buildup  of  the  angle  of  drift  and  the 
result  is  an  eventual  stable  limit  cycle. 

The  ocean  environment  can  be  expected  to  introduce  many  combinations 
of  disturbance  forces,  which  may  or  may  not  be  periodic.  A  preferred  method 
for  simulating  many  disturbances  such  as  sensor  noise,  ocean  current 
fluctuation,  and  vehicle  acceleration  fluctuations  is  to  model  them  as  'white 
noise'.  Figure  29  is  presented  to  demonstrate  the  effects  of  including 
constant,  zero-mean  disturbances  in  sway  and  yaw  on  the  marginally  stable 
system  of  configuration  'B '.  The  disturbances  are  developed  using  MATLAB's 
random  number  generator  with  a  uniform  distribution.  Che  rnndom  numbers 
are  then  scaled  to  simulate  values  that  may  be  expected.  The  variance  of  the 
sway  and  yaw  acceleration  disturbances,  respectively,  for  Figure  29  are  0.003 
(ft/s2)2  and  0.002  (rad/s2)2.  The  resulting  simulation  bears  little  resemblance 
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to  Figure  24,  although  the  initial  conditions  and  values  for  Xg  and  Zg  are  the 
same.  With  even  these  relatively  minor  disturbances  acting  on  the  system  of 
configuration  'B\  a  rather  large,  non-symmetric  oscillation  in  both  roll  and 
angle  of  drift  is  evident.  Although  the  system  is  still  stable,  with  the  mean  of 
both  the  roll  angle  and  angle  of  drift  equalling  zero,  a  limit  cycle  similar  to 
that  of  the  unstable  configuration  (Figure  25)  has  developed.  As  the  angle  of 
drift  fluctuates  between  positive  three  degrees  and  negative  four  degrees,  the 
angle  of  roll  varies  between  positive  and  negative  two  degrees.  Increasing 
the  scaling  (which  increases  the  variance)  for  the  disturbances  would  be 
expected  to  increase  the  fluctuations  until  stability  is  lost.  Similarly,  it  is 
true  that  the  small  disturbances  acting  in  Figure  29  have  a  greatly  reduced 
effect  on  the  system  of  configuration  'A',  which  has  greater  stability. 
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IV.  PERTURBATION  ANALYSIS 

A.  BACKGROUND 

The  previous  section  demonstrates  that  knowledge  of  the  separation 
between  longitudinal  centers  of  buoyancy  and  gravity  is  critical  in 
determining  system  stability.  If  this  quantity  is  the  dominant  factor  in  a 
stability  analysis,  then  an  approximation  for  the  degree  of  stability  may  be 
developed  by  applying  perturbation  theory.  Although  perturbation  results 
are  only  approximations,  their  advantage  over  numerical  methods  is  in 
illustrating  the  degree  to  which  a  solution  depends  on  the  variable(s) 
involved.  From  the  fundamental  theorem  of  perturbation  theory  [Ref.  6],  we 
seek  a  solution  to  the  characteristic  equation  AX4  +  BX  +  CX  +  DX  +  E  =  0  of 
the  form: 

n=oo 

Me)  =    X  an  en  (a  i) 

n=0  ' 

where  e  is  the  variable  of  interest  and  the  solution  approximates  the 
numerical  solution  in  the  region  where  e  is  smaEl.  The  constant  coefficients 
(ao,  ai,  ...,  an)  axe  all  independent  of  £,  and  are  alltequal  to  zero  for  small  e. 

B.  PERTURBATION  FORMULATION 

Recalling  Equations  2.4f  -  2.4j,  which  defhned  the  coefficients  of  the 
quartic  characteristic  equation,  it  would  obviously  be  desirable  to  simplify  the 
equations  as  much  as  possible  prior  to  formulating  the  perturbation 
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approximation.        The     variable     of    interest    will     be    defined     as: 

e  =  Xg  -  Xb  (4.2) 

where  the  nominal  value  for  Xb  is  zero  and  the  perturbation  will  be 

performed  around  e  =  0.  By  comparing  Figures  30  and  31  it  is  clear  that  by 
allowing  the  hydrodynamic  coefficients  Kv,  Kf-,  Yp,  Np,  Yp,  and  1^>  to    equal 

zero  a  very  good  approximate  solution  to  Equation  (2.4e)  is  obtained  for  both 
configurations  'A'  and  'B\  This  simplification  reduces  the  descriptor 
coefficients  T ,  T,  'o',  and  'p'  from  Table  1  to  zero,  which  simplifies  the 
coefficients  of  the  coupled  equations  of  motion  to  Equations  (4.3a)  -  (4.3e) 
below. 

A  =  (Ixx-  Kp)[(M-Yv)(Izz-N,)-(MXg-Y,)(MXg-Nv)]  + 

(MZg)2(Izz-Nr)  (4.3a) 

B  =  (KvU)(MZgXIzz-N,)  +  (MZg)2(Nr-MXg)(U)- 

(MZg)(MZg  +  Kr)(MXg-Nv)(U)-axx-Kp)[(M-Yv)(Nr-MXg)(U)  + 
(YvUXIzz-Nr)-(NvU)(MXg-Yr)-(MXg-Nv)(Yr-M)(U)]- 
(KpU)[(M-Yv)(Izz-Nr)-(MXg-Yr)(MXg-Nv)]  (4.3b) 

C  =  (Ixx-Kp)[(YvU2XNr-MXg)-(NvU2)(Yr-M)]  + 
(ZgWX(M- Yv)(Izz-Nr)-(MXg  -Yr)(MXg  -Nv)]  + 
(KpU)[(M-Yv)(Nr-MXg)(U)  +  (YvU)(I„  -Nr)  - 
(NyUXMXg-Yt)  -  (MXg-Nv)(Yr-M)(U)]  + 
(MZg)[(XgV/XMXg  -Y,  WNvU2)(MZg  +  Kr)  +  (KvU2XNr  -MXg)] 


D  =    -(ZgW)[(M-Yv)(Nr-MXg)(U)  +  (YvU)(Izz-N,)- 
(NvUXMXg  -  Yr)  -  (MXg  -  Nv)(Yr  - M)(U)]  - 
(KpU3)[(Yv)(Nr-MXg)-(Nv)(Yr-M)]  +  (KvU)(XgW)(MXg-  Yr)- 
(MZg)(XgW)(U)(Yr  -  M)  -  (MZg  +Kr)(U)(XgW)(M  -Yv)  (4  3d) 

E  =  (ZgW)[(YvU2)(Nr-MXg)-(NvU2)(Yr-M)]  -  (KvXgWU2)(Yr -M)  + 

(YvXgWU2)(MZg+Kr)  (4>3e) 


Recalling  the  definitions  of  the  coefficients  for  the  uncoupled  systems: 

AL   =  (M-Yv)(Izz-Nr)  -  (MXg-Yr)(MXg-Nv) 

BL  =    -[(Izz-Nr)(YvU)  +  (M-Yv)(Nr-MXg)(U)  + 
(Yr  -M)(Nv  -MXg)(U)  +  (Yr  -  MXg)(NvU)] 

CL  =  (YvU2XNr-MXg)  -  (NvU2XYr-M) 

AR  =  Ixx-Kp  BR  =  -KpU  CR  =  ZgW 

Substituting  in  Equations  (4.3a)  -  (4.3e)  yields: 

A  =  ARAL  +  (MZg)2(Izz-Nr)  (4.4a) 

B  =  ARBL  +  BRAL  +  (MZg)2(Nr-MXg)(U)  +  (4>4b) 

(MZgXdzz  -Nr)(KvU)-(MZg  +Kr)(MXg  -NV)(U)] 


C  =  ArCl  +  CrAl  +  BrBl  +  (MZg)[(XgW)(MXg-Yr)  + 
(KvU2XNr  -  MXg)  -  (NvU2)(MZg  +  Kr)] 


(4.4c) 


D  =  CrBl  +  CLBR  +  (KvU)(XgW)(MXg-Yr)  - 

(MZg  +  KrXXgW)(U)(M-Yv)  -  (MZg)(XgW)(U)(Yr-M)  (4.4d) 

E  =  CrCl  +  (YvU2)(XgWXMZg  +  Kr)  -  (KvXgWU2)(Yr -M)  (4  4g) 
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In  order  to  further  simplify  these  coefficients,  terms  of  order  (Zg)    and  (Xg) 

will  be  neglected;  the  effects  of  doing  this  are  small  and  may  be  seen  in 

Figures  32  and  33.  This  reduces  the  coefficients  to: 

A  =  ARAL  (4.5a) 

B  =  ARBL  +  BRAL  +  a  (Xg)  +  Kl  (4.5b) 

C  =  ArCl  +  CrAl  +  BRBL  +  (3  (Xg)  +  K2  (4.5c) 

D  =  CRBL  +  CLBR  +  Y(Xg)  (4.5d) 

E  =  CRCL  +  6(Xg)  (4.5e) 

where:  a  =  -M2ZgKrU  (4.6a) 

(3  =    -(MZg)[(WYr)  +  (MKvU2)]  (4.6b) 

Y  =  (WUXMZg(Yv-Yr)+  KrYv-KrM-KvYfl  (4.6c) 

8  =  (WU^tCYyXMZg+Kr)  -  (Kv)(YrM)]  (4.6e) 

Kl  =  (MZg)[(Izz-Nr)(KvU)  +  (NvKrU)]  (4.6£) 

K2  =  (MZgU2)(NrKv  -  NvK,-)  (4.6g) 

Carrying  through  the  computations  results  in  the  following  expressions: 

(ArAl)  ao4  +  (ArBl+BrAl+K1)  a03  +  (ArCl+CrAl+BrBl+K2)  a02  + 

(CrBl+ClBr)  a0  +  CRCL  =  0  (4.7a) 

_  -(aa03  +  (3a02  +  ya0  +  8) __  (4  ?1. 

1        (4A)a03+3(B-aXg)ao2+2(C-pXg)a0  +  (D-YXg) 

-ai2(C-pX8)  -  ai2a0(6A+B-aXg)  -  Klapai2  -  3oeao2ai-2fiaoai-Yai 
&2  "  (4A)a03  +  3(B-aXg)ao2  +  2(C-|3Xg)ao+(D-YXg) 

(4.7c) 
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The  characteristic  equation  AX4  +  BX  +  CX  +  DX  +  E  =  0  may  now  be 
expressed  in  terms  of  a  second  order  perturbation  which  depends  on  a  single 
variable  e  as: 

X(e)  =  ao  +  ai  e  +  a2  e2  (4.8) 

Figures  34  through  37  depict  the  results  for  both  first  and  second  order 
perturbation  analysis  about  Xg  close  to  zero.  Results  for  both  configurations 

'A'  and  'B'  demonstrate  that  this  method  for  representing  the  degree  of 
stability  in  the  vicinity  of  a  nominal  (Xg  -  Xb)  is  quite  satisfactory,  provided 

the  roots  have  no  complex  conjugates.  In  the  region  where  the  roots  become 
complex,  the  problem  is  no  longer  a  regular  perturbation,  and  methods  to 
solve  the  singular  perturbation  must  be  developed. 
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V.    CONCLUSIONS  AND  RECOMMENDATIONS 

When  considering  an  inherently  stable  situation  for  an  analysis  of  sway, 
yaw,  and  roll  response,  the  linear  simulations  compare  favorably  with 
the  nonlinear  simulations.  For  a  marginally  stable  or  unstable  system, 
only  the  nonlinear  simulation  can  predict  the  existence  and  extent  of  any 
non-zero  steady  states  or  limit  cycles  that  may  occur. 

Obviously,  in  the  design  of  a  submersible,  careful  attention  must  be  paid 
to  selection  of  both  the  metacentric  height  and  the  separation  between  the 
longitudinal  centers  of  gravity  and  buoyancy.  Instability  associated  with 
the  dynamic  coupling  effects  can  be  minimized  by  increasing  the 
metacentric  height.  The  uncoupled  system  stability  predictions  are  not 
reliable  when  there  is  separation  between  Xg  and  Xb. 

It  can  be  concluded  from  the  analysis  that  the  coupling  of  roll  into  sway 
and  yaw  for  the  linearized  equations  of  motion  is  not  very  significant 
when  Xg  equals  Xb. 

Recommendations  for  future  modelling  research  include  an  expansion 
to  incorporate  coupling  effects  between  the  vertical  and  horizontal  planes. 
It  is  well  known  that  the  coupling  between  translation  and  attitude  is  a 
severe  limitation  in  the  design  of  submersibles. 

The  effects  of  varying  hydrodynamic  derivatives  and  initial  conditions,  as 
well  as  incorporating  disturbances  in  the  nonlinear  simulations  deserve 
more  attention.  Additionally,  the  complexity  of  the  nonlinear  system 
resulting  from  including  the  vertical  plane  is  greatly  increased.  More 
studies  of  the  nonlinear  system  resulting  from  this  union  are  also 
required. 
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APPENDIX  A 

NONLINEAR  SIMULATION  PROGRAM: 

**************************************************** 

%  THIS  PROGRAM  USES  HAMMING'S  METHOD  TO  SOLVE  A  SYSTEM 
%  OF  NONLINEAR  EQUATIONS.  SIMILAR  TO  A  FOURTH  ORDER 
%  RUNGE-KUTTA  IN  TERMS  OF  LOCAL  ERROR,  IT  REQUIRES  ONLY 
%  TWO  FUNCTION  CALCULATIONS  PER  STEP  VICE  FOUR. 

function[v,phi]  =  nonlinsim_hamming(UO,Xg,Zg) 

%  VEHICLE  PARAMETERS  FOLLOW 

W  =  12000.0;  Ixx  =  1760.0;  Iyy  =  9450.0;  Izz=10700.0;  L  =  17.425; 
RHO=1.94;  G=32.2;  M  =  W/G;  BUOY  =  W;  Zb=0.0;   Xb  =  0.0; 
ZZ  =  0.5*RHO*LA2; 

%  SWAY  HYDRODYNAMIC  COEFFICIENTS 


Ypdot 

=  1.270E-04*ZZ*LA2; 

Yrdot 

=  1.240E-03*ZZ*LA2; 

Ypq 

=  4.125E-03*ZZ*LA2; 

Yqr 

=  -6.510E-03*ZZ*LA2; 

Yvdot 

=  -5.550E-02*ZZ*L; 

Yp 

=  3.055E-03*ZZ*L; 

Yvq 

=  2.360E-02*ZZ*L; 

Ywp 

=  2.350E-01*ZZ*L; 

Ywr 

=  -1.880E-02*ZZ*L; 

Yv 

=  -9.310E-02*ZZ; 

Yvw 

=  6.840E-02*ZZ; 

%  NOMINAL  VALUE  FOR  Yr  =  2.970E-02*ZZ*L 
%  CONFIGURATION  'A':  Yr  =  -3.500E-02*ZZ*L 
%  CONFIGURATION  'B':    Yr  =  -5.940E-02*ZZ*L 

CC  =  inputCenter  model  configuration  choice;  either  1  for  A  or  2  for  B') 


if  CC<2 

Yr  =  -3.500e-02*ZZ*L; 
else 

Yr  =  -5.940e-02*ZZ*L; 
end 

ROLL  HYDRODYNAMIC  COEFFICIENTS 
Kpdot   =  -1.01E-03*ZZ*LA3;               Krdot    = -3.37E-05*ZZ*LA3; 
Kpq       = -6.93E-05*ZZ*LA3;               Kqr       =  1.68E-02*ZZ*LA3; 
Kvdot    =  1.27E-04*ZZ*LA2;               Kp        = -1.10E-02*ZZ*LA2; 
Kr         = -8.41E-04*ZZ*LA2;                Kvq       = -5.115E-03*ZZ*LA2; 
Kwp      =  -1.27E-04*ZZ*LA2;                Kwr      =  1.39E-02*ZZ*LA2; 
Kv         =  3.055E-03*ZZ*L;                  Kvw      =  -1.87E-01*ZZ*L; 
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YAW  HYDRODYNAMIC  COEFFICIENTS 


Npdot 

=  -3.370E-05*ZZ*LA3; 

Nrdot 

=  -3.400E-03*ZZ*LA3; 

Npq 

=  -2.110E-02*ZZ*LA3; 

Nqr 

=  2.750E-03*ZZ*LA3; 

Nvdot 

=  1.240E-03*ZZ*LA2; 

Np 

=  -8.405E-04*ZZ*LA2; 

Nr 

=  -1.640E-02*ZZ*LA2; 

Nvq 

=  -9.990E-03*ZZ*LA2; 

Nwp 

=  -1.750E-02*ZZ*LA2; 

Nwr 

=  7.350E-03*ZZ*LA2; 

Nv 

=  -1.484E-02*ZZ*L; 

Nvw 

=  -2.670E-02*ZZ*L; 

%  THE  FOLLOWING  ARE  USED  FOR  EVALUATING  THE  INTEGRALS 
%  IN  THE  SWAY  AND  YAW  EQUATIONS  OF  MOTION.   <X'  IS  THE 
%  DISTANCE  IN  FEET  ALONG  THE  LONGITUDINAL  AXIS  AND  'HH' 
%  IS  THE  VEHICLE  HEIGHT.  ALL  VALUES  ARE  IN  FEET. 


X(l)=-105.9/12; 

X(4)=-94.3/12; 

X(7)=-66.3/12; 

X(10)=79.2/12; 

X(13)=91.2/12; 

X(16)=101.2/12; 


X(2)=-104.3/12; 

X(5)=-87.3/12; 

X(8)=-55.8/12; 

X(ll)=83.2/12; 

X(14)=95.2/12; 

X(17)=102.1/12; 


X(3)=-99.3/12; 

X(6)=-76.8/12; 

X(9)=72.7/12; 

X(12)=87.2/12; 

X(15)=99.2/12; 

X(18)=103.2/12; 


HH(1)=0.0/12; 

HH(4)=13.96/12; 

HH(7)=29.36/12; 

HH(10)=30.00/12; 

HH(13)=21.44/12; 

HH(16)=9.12/12; 


HH(2)=2.28/12; 

HH(5)=19.76/12; 

HH(8)=31.85/12; 

HH(11)=27.84/12; 

HH(14)=17.12/12; 

HH(17)=6.72/12; 


HH(3)=8.24/12; 

HH(6)=25.1/12; 

HH(9)=31.85/12; 

HH(12)=25.12/12; 

HH(15)=12.0/12; 

HH(18)=0.00/12; 


%  THE  INITIAL  CONDITIONS  ARE  SET  PRIOR  TO  INTEGRATION 
%  BY  CALLING  A  LINEAR  INTEGRATION  PROGRAM  NAMED 
%  'EULER'.  THE  ERRORS  ASSOCIATED  WITH  USING  LINEAR 
%  SOLUTIONS  FOR  THE  FIRST  FIVE  TIME  INTERVAL  STEPS  ARE 
%  SMALL  AND  DO  NOT  AFFECT  THE  NONLINEAR  SOLUTIONS. 

[p,pdot,v,vdot,r,rdot,phi,phidot]  =  euler(dt); 

p(l:5)=p; 

pdot(l:5)=pdot; 

pdotmod(5)=pdot(5); 

v(l:5)=v; 

vdot(l:5)=vdot; 

vdotmod(5)=vdot(5); 
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Kl:5)=r; 

rdot(l:5)=rdot; 

rdotmod(5)=rdot(5); 

phi(l:5)=  phi; 
phidot(l:5)=  phidot; 
phidotmod(5)=phidot(5); 

%  THE  TIME  INTERVAL  IS  'dt\  THE  FINAL  TIME  IS  'tfinal',  AND  A 
%  VALUE  CALLED  'stopnumber'  IS  SET  TO  ALLOW  ACCESS  TO  THE 
%  DATA  AFTER  THIS  NUMBER  OF  ITERATIONS.  THE  VALUE  FOR 
%  'stopnumber  MAY  BE  CHANGED  WHILE  THE  KEYBOARD  IS 
%  ACTIVE  TO  ALLOW  SUBSEQUENT  PROGRAM  INTERACTION. 
%  RETURN  TO  THE  PROGRAM  IS  ACHIEVED  BY  TYPING  'return' 
%  FOLLOWED  BY  THE  'enter'  KEY. 

dt  =  inputCenter  the  time  interval  step  size') 

tfinal  =  inputCenter  the  final  time') 

stopnumber  =  inputCenter  the  value  for  stopnumber') 

%  THIS  SECTION  ALLOWS  FOR  RANDOM  DISTURBANCES  IN  SWAY 
%  AND  YAW 

rand('uniform') 

vdotdist  =  rand(l,(tfinal/dt)+l); 

vdotdist  =  0.05*(vdotdist  -  mean( vdotdist)); 

rdotdist  =  rand(l,(tfinal/dt)+l); 

rdotdist  =  0.04*(rdotdist  -  mean(rdotdist)); 

%  THIS  SECTION  PROVIDES  FOR  CONTINUATION  OF  SOLUTIONS. 
%  IF  THE  MATRICES  BECOME  VERY  LARGE,  IT  IS  RECOMMENDED 
%  THAT  THE  CURRENT  VALUES  BE  SAVED  AND  A  NEW 
%  SIMULATION  BEGUN  AS  MATLAB  SLOWS  NOTICEABLY  WITH 
%  LARGE  MATRICES. 
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vcorr(5)=  ;  rcorr(5)= 

pcorr(5)=  ;  phicorr(5)= 


%    vpred(5)=  ;  rpred(5)= 

%     ppred(5)=  ;  phipred(5)= 

%     pdotmod(5)=     ;  rdotmod(5)= 
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%   COMMENCE  INTEGRATION 

forj  =  6:(tfinaI/dt)  +  l; 

%  THIS  SECTION  USES  A  HAMMING  PREDICTOR-CORRECTOR  TO 
OBTAIN  VALUES 

vpred(j)=v(j-4)+(4*dt/3)*(2*vdot(j-l)-vdot(j-2)+2*vdot(j-3)); 
rpred(j)=r(j-4)+(4*dt/3)*(2*rdot(j-l)-rdot(j-2)+2*rdot(j-3)); 
ppred(j)=p(j-4)+(4*dt/3)*(2*pdot(j-l)-pdot(j-2)+2*pdot(j-3)); 
phipred(j)=phi(j-4)+(4*dt/3)*(2*p(j-l)-p(j-2)+2*p(j-3)); 


ifj<7 

vmod(j)  = 
rmod(j)  = 
pmod(j)  = 
phimod(j): 

else 

vmod(j)= 
rmod(j)= 
pmod(j)= 
phimod(j)= 

end 


vpred(j); 
rpred(j); 
ppred(j); 
phipred(j); 

vpred(j  )-(l  12/12  l)*(vpred(j-l)-vcorr(j-l)); 
rpred(j)-(112/121)*(rpred(j-l)-rcorr(j-l)); 
ppred(j)-(112/121)*(ppred(j-l)-pcorr(j-l)); 
phipredCJ )-( 1 12/12 1  )*(phipred(j- 1  )-phicorr(j- 1 )); 


tmplmod=vmod(j);  tmp2mod=rmod(j); 
[SWAYMOD,YAWMOD]=crossflowintmod(tmplmod,tmp2mod,X,HH); 

vdotmod(j)=        (l/(M-Yvdot))*((Ypdot+M*Zg)*pdotmod(j-l)+... 
Yp*UO*pmod(j)+(Yrdot-M*Xg)*rdotmod(j-l)+... 
Yv*U0*vmod(j)+(Yr-M)*U0*rmod(j)+SWAYMOD); 

vcorr(j)=  0.125*(9*v(j-l)-v(j-3)+3*dt*(vdotmod(j)+2*vdot(j-l)-... 

vdot(j-2))); 

rdotmod(j)=        (l/(Izz-Nrdot))*(Npdot*pdotmod(j-l)+... 

(Nvdot-M*Xg)*vdotmod(j  )+Np*U0*pmod(j )+. . . 

(Nr-M*Xg)*UO*rmod(j)+Nv*UO*vmod(j)+... 

Xg*W*sin(phimod(j))+YAWMOD); 

rcorr(j)=  0.125*(9*r(j-l)-r(j-3)+3*dt*(rdotmod(j)+2*rdot(j-l)-... 

rdot(j-2))); 

pdotmod(j)=        (l/(Ixx-Kpdot))*((Kvdot+M*Zg)*vdotmod(j)+... 
Krdot*rdotmod(j)+Kp*UO*pmod(jW... 
(Kr+M*Zg)*UO*rmod(j)+Kv*UO*vmod(j)-... 
Zg*W*sin(phimod(j))); 
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pcorr(j)=  0.125*(9*p(j-l)-p(j-3)+3*dt*(pdotmod(j)+2*pdot(j-l)-... 

pdot(j-2))); 

phicorr(j)=  0.125*(9*phi(j-l)-phi(j-3)+3*dt*(pmod(j)+2*p(j-l)-... 

p(j-2))); 

v(j)=  vcorr(j)+(9/121)*(vpred(j)-vcorr(j)); 

r(j)=  rcorr(j)+(9/121)*(rpred(j)-rcorr(j)); 

p(j)=  pcorr(j)+(9/121)*(ppred(j)-pcorr(j)); 

phi(j)=  phicorr(j)+(9/121)*(phipred(j)-phicorr(j)); 

beta(j)  =  -atan(v(j)/5)*180/pi; 

%  THIS  SECTION  HAS  THE  NON-LINEAR  EQUATIONS  OF  MOTION. 
%  REMEMBER  TO  REMOVE  THE  DISTURBANCES  FROM  THE  SWAY 
%  AND  YAW  DERIVATIVES  IF  THEY  ARE  NOT  DESIRED. 

tmpl=v(j);  tmp2=r(j); 
[SWAY.YAW]  =  crossflowint(tmpl,tmp2,X,HH); 

vdot(j)=  ((l/(M-Yvdot))*((Ypdot+M*Zg)*pdot(j-l)+Yp*UO*p(j)  +... 
(Yrdot-M*Xg)*rdot(j-l)+Yv*UO*v(j)+(Yr-M)*UO*r(j)+... 
SWAY))+vdotdist(j); 

rdot(j)=     ((l/(Izz-Nrdot))*(Npdot*pdot(j-l)+(Nvdot-M*Xg)*vdot(j)+... 
Np*UO*p(j)+(Nr-M*Xg)*UO*r(j)+Nv*UO*v(j)+... 
Xg*W*sin(phi(j))+YAW))+rdotdist(j); 

pdot(j)=    (l/(Ixx-Kpdot))*((Kvdot+M*Zg)*vdot(j)+Krdot*rdot(j)+... 
Kp*UO*p(j)+(Kr+M*Zg)*UO*r(j)+Kv*UO*v(j)-... 
Zg*W*sin(phi(j))); 

if  j  ==  stopnumber 
keyboard 
end 

end 
return 
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LINEAR  SIMULATION  PROGRAM: 

************************************************************************ 

%  THIS  PROGRAM  IS  THE  LINEAR  EQUATION  OF  MOTION 

%  INTEGRATION  EMPLOYING  A  SIMPLE  EULER  METHOD  TO 

%  OBTAIN  STARTING  VALUES  FOR  THE  NONLINEAR  SIMULATION. 

function[p,pdot,v,vdot,r,rdot,phi,phidot]  =  euler(dt) 

%  THE  MATRIX  'A  CORRESPONDS  TO  THE  DAMPING  MATRIX  AND 
%  MATRIX  'B'  CORRESPONDS  TO  THE  MASS  MATRIX. 


A(l,l)  =  Kp*U0 
A(l,3)  =  Kv*U0: 
A(2,l)=1.0; 
A(2,3)  =  0.0; 
A(3,l)  =  Yp*U0 
A(3,3)  =  Yv*U0 
A(4,l)  =  Np*U0 
A(4,3)  =  Nv*U0 


A(l,2)  =  -((Zg*WHZb*B)); 
A(l,4)  =  (Kr+(M*Zg))*U0; 
A(2,2)  =  0.0; 
A(2,4)  =  0.0; 
A(3,2)  =  0.0; 
A(3,4)  =  (Yr-M)*U0; 
A(4,2)  =  ((Xg*W)-(Xb*B)); 
A(4,4)  =  (Nr-(M*Xg))*U0; 


B(l,l)  =  Ixx-Kpdot; 
B(l,3)  =  -(Kvdot+(M*Zg)); 
B(2,l)  =  0.0; 
B(2,3)  =  0.0; 

B(3,l)  =  -(Ypdot+(M*Zg)); 
B(3,3)  =  M-Yvdot; 
B(4,l)  =  -Npdot; 
B(4,3)  =  (M*Xg)-Nvdot; 

C  =  [inv(B)*A]; 

%  INITIAL  CONDITIONS 


B(l,2)  =  0.0; 

B(l,4)  =  -Krdot; 

B(2,2)  =  1.0; 

B(2,4)  =  0.0; 

B(3,2)  =  0.0; 

B(3,4)  =  (M*Xg)-Yrdot; 

B(4,2)  =  0.0; 

B(4,4)  =  Izz-Nrdot; 


timed)  =  0.0; 
p(l)  =  0.0; 
phi(l)  =  1.0*pi/180; 
psi(l)  =  0.0; 
v(l)  =  0.0; 
r(l)  =  0.0; 
x(l)  =  0.0; 
xdot(l)  =  0.0; 
tfinal  =  1; 


pdotd)  =  0.0; 
phidot(l)  =  0.0; 
psidot(l)  =  0.0; 
vdot(l)  =  0.0; 
rdot(l)  =  0.0; 
y(l)  =  0.0; 
ydot(l)  =  0.0; 
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%  THIS  SECTION  ALLOWS  FOR  RANDOM  DISTURBANCES  IN  SWAY 
%  AND  YAW  FOR  THE  LINEAR  EQUATIONS  AS  WELL.  MATLAB 
%  ALLOWS  FOR  EITHER  A  UNIFORM  OR  NORMAL  DISTRIBUTION 
%  OF  RANDOM  NUMBER. 

rand(  uniform') 

vdotdist  =  rand(l,(tfinal/deltat)+l); 
vdotdist  =  0.05*(vdotdist  -  mean(vdotdist)); 
rdotdist  =  rand(l,(tfinal/deltat)+l); 
rdotdist  =  0.04*(rdotdist  -  mean(rdotdist)); 

%  COMMENCE  EULER  INTEGRATION 

for  j  =  2:(tfinal/deltat)  +  1; 

P(j)  =  P(j-D  +  pdot(j-l)*deltat; 

phi(j)=  phi(j-l)  +  phidot(j-l)*deltat; 

psi(j)=  psi(j-l)  +  psidot(j-l)*deltat; 

v(j)  =  v(j-l)  +  vdot(j-l)*deltat; 

beta(j)=  -atan(v(j)/U0); 

r(j)  =  r(j-l)  +  rdot(j-l)*deltat; 

x(j)=  x(j-l)  +  xdot(j-l)*deltat; 

y(j)=  y(j-l)  +  ydot(j-l)*deltat; 


pdot(j)  = 

phidot(j)= 

vdot(j)  = 

rdot(j)  = 

psidot(j) : 

xdot(j)  = 

ydot(j)  = 

time(j)  = 

end 
return 


C(l,l)*p(j)  +  C(l,2)*phi(j)  +  C(l,3)*v(j)  +  C(l,4)*r(i) 

C(2,l)*p(j)  +  C(2,2)*phi(j)  +  C(2,3)*v(j)  +  C(2,4)*r(j) 

C(3,l)*p(j)  +  C(3,2)*phi(j)  +  C(3,3)*v(j)  +  C(3,4)*r(j) 

C(4,l)*p(j)  +  C(4,2)*phi(j)  +  C(4,3)*v(j)  +  C(4,4)*r(j) 

r(j)*cos(phi(j)); 

U0*cos(psi(j))  -  v(j)*sin(psi(j))*cos(phi(j)); 

UO*sin(psi(j))  +  v(j)*cos(psi(j))*cos(phi(j)); 

(j*deltat)  -  deltat; 
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CROSSFLOW  INTEGRAL  PROGRAM: 

************************************************************************ 

%  THIS  IS  A  NUMERICAL  INTEGRATION  PROGRAM  TO  CALCULATE 
%  THE  DRAG  FORCES  IN  THE  HORIZONTAL  PLANE  UTILIZING  THE 
%  TRAPEZOIDAL  RULE.   THE  CALL  TO  'crossflowintmod'  IS 
%  IDENTICAL,  ONLY  USING  DIFFERENT  VARIABLE  NAMES  AND 
%  VALUES  FOR  THE  MODIFICATION  PORTION  OF  HAMMING'S 
%  METHOD. 

function[SWAY,YAW]=crossflowint(tmpl,tmp2,X,HH) 

RHO  =  1.94;     CDy  =  0.35;     SWAY  =  0.0;     YAW  =  0.0; 
v  =  tmpl;  r  =  tmp2; 

fork  =1:18; 
ifabs(v+X(k)*r)<le-6 
UCF(k)  =  0.0; 
else 

UCF(k)    =  (v+X(k)*r)/(abs(v+X(k)*r)); 
end 
CFLOW(k)  =  CDy*HH(k)*((v+X(k)*r)A2); 
SWAYl(k)  =  CFLOW(k)*UCF(k); 
YAWl(k)  =  CFLOW(k)*X(k)*UCF(k); 
end 

jj  =  1:17; 

SWAY  =  -0.25*RHO*sum((SWAYl(jj)+SWAYl(jj+l)).*(X(jj+l)-X(jj))); 

YAW  =  -0.25*RHO*sum((YAWl(jj)+YAWl(jj+l)).*(X(jj+l)-X(jj))); 
return 


EIGENVALUE  CALCULATION  PROGRAM; 

************************************************************************ 

%  THIS  PROGRAM  DETERMINES  THE  FOUR  ROOTS  FOR  THE  ROLL 
%  COUPLED  EQUATIONS  OF  MOTION  OVER  A  RANGE  OF  Xg. 

function[all_roots,xg]  =  findroots(UO,Xgmin,stepXg,Xgmax) 


W=  12000.0;    Ixx  =  1760.0;    Iyy  =  9450.0;    Izz=10700.0; 
L=  17.425;       RHO=1.94;       G=32.2;  M  =  W/G; 

BUOY  =  W;Zb=0.0; 
ZZ  =  0.5*RHO*LA2; 
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%  SWAY  HYDRODYNAMIC  COEFFICIENTS 

Ypdot=  1.27E-04*ZZ*LA2;  Yrdot=  1.24E-03*ZZ*LA2; 
Yvdot  =  -5.55E-02*ZZ*L;  Yp      =  3.055E-03*ZZ*L; 

Yv      =  -9.31E-02*ZZ; 

CC  =  inputCenter  model  configuration  choice;  either  1  for  A  or  2  for  B') 

if  CC<2 

Yr  =  -3.500e-02*ZZ*L; 
else 

Yr  =  -5.940e-02*ZZ*L; 
end 

%  ROLL  HYDRODYNAMIC  COEFFICIENTS 

Kpdot  =  -1.01E-03*ZZ*LA3;  Krdot  =  -3.37E-05*ZZ*LA3; 
Kvdot  =  1.27E-04*ZZ*LA2;  Kp  =  -1.10E-02*ZZ*LA2; 
Kr        =  -8.41E-04*ZZ*LA2;     Kv       =  3.055E-03*ZZ*L; 

%  YAW  HYDRODYNAMIC  COEFFICIENTS 

Npdot  =  -3.370E-05*ZZ*LA3;  Nrdot  =  -3.400E-03*ZZ*LA3; 
Nvdot  =  1.240E-03*ZZ*LA2;  Np  =  -8.405E-04*ZZ*LA2; 
Nr     =  -1.640E-02*ZZ*LA2;      Nv       =  -1.484E-02*ZZ*L; 

rows=(abs(Xgmax-Xgmin)/stepXg)+  1; 
all_roots  =  zeros(rows,4);        xg  =  zeros(rows,l); 

Zg  =  inputCenter  value  for  Zg') 

a  =  (Ixx-Kpdot);  b  =  (Kp*U0);  c  =  (Zg*W)-(Zb*BUOY); 

d  =  (M*Zg)+Kvdot;  e  =  Kv*U0;  f  =  Krdot; 

g  =  (M*Zg*U0)  +  (Kr*U0);  h  =  Ypdot  +  (M*Zg); 

i  =  Yp*U0;  j  =  M  -  Yvdot;  k  =  Yv*U0; 

0  =  Npdot;  u  =  Izz  -  Nrdot;        p  =  Np*U0; 
m  =  (Yr*U0)  -  (M*U0);  x  =  Nv*U0; 

Xg  =  Xgmin:stepXg:Xgmax; 
Xb  =  inputCenter  either  0.0  or  Xg  for  value  of  Xb); 

1  =  (M*Xg)  -  Yrdot;  q  =  (Xb.*BUOY)-(Xg*W); 
r  =  (M*Xg)-Nvdot;  w  =  (Nr*U0)-(M*Xg*U0); 

A  =  (((a).*(j+u-l.*r)+(d).*(u*h+o.*l)-(f).*(r.*h-rO*j));; 

B  =  ((e).*(u*h+o.*l)-(d).*(u*i-w.*h+o.*m-p.*l)-(a).*... 
(j.*w+k*u-l.*x-r.*m)-(b).*(j*u-l.*r)+(f).*(r.*i-x... 
.*h+o*k-p*j)-(g).*(r.*h+o*j)); 
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C  =  ((a).*(k.*w-m  *x)+(b)  *(j.*w+k*u-l.*x-r.*m)+(c).*... 
(j*u-l.*r)-(d).*(q.*l-w.*i+m.*p)-(e).*(u*i-w.*h+... 
o.*m-p.*l)+(f).*(q.*j-x*i+p*k)+(g).*(r.*i-x.*h+o*k... 

-P.*j»; 

D  =  ((g).*(q.*j-x.*i+p*kHf.*q*k>+(d  *q.*m>-(e).*(q.*l... 

-w.*i+m.*pHc).*(j.*w+k*u-l.*x-r*m)-(b).*(k*w-m.*x)); 

E  =  ((c)  *(k  *w-m.*x)+(e.*q.*m)-(g.*q.*k)); 

z  =  [A*  B'  C  J)'  El; 

j  =  lrrows;  poly(j,:)=z(j,:);  rts  =  roots(poly(jj,:))'; 

for  jj  =  l:rows; 
all_roots(jj,l:4)  =  eval(rts)'; 
end 

xg=Xg'; 

end; 

return 


FIRST  AND  SECOND  ORDER  PERTURBATION 
APPROXIMATION  COMPUTER  PROGRAM, 

***************************************************************** 

%  THIS  PROGRAM  DEVELOPS  THE  1ST  &  2ND  ORDER  PERTURBATION 
%  APPROXIMATION  SOLUTIONS  FROM  SECTION  IV 

function[xlocation,perturbrootl,perturbroot2,xg]  =perturbanalysis(UO) 

W  =  12000.0;         Ixx  =  1760.0;  Iyy  =  9450.0;  Izz=10700.0; 

L=  17.425;  RHO=1.94;  G=32.2;  M  =  W/G; 

BUOY  =  W;  Zb=0.0;  Xb  =  0.0; 

ZZ  =  0.5*RHO*LA2; 

%  SWAY  HYDHODYNAMIC  COEFFICIENTS 

Ypdot=  0.0;  Yrdot^  1.24E-03*ZZ*LA2; 

Yvdot  =  -5.55E-02*ZZ*L;  Yp    =0.0; 

Yv    =  -9.31E-02*ZZ; 

CC  =  inputCenter  model  configuration  choice;  either  1  for  A  or  2  for  B') 

if  CC<2 

Yr  =  -3.500e-02*ZZ*L: 
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else 

Yr  =  -5.940e-02*ZZ*L; 

end 

%  ROLL  HYDRODYNAMIC  COEFFICIENTS 
Kpdot  =  -1.01E-03*ZZ*LA3;  Krdot  =  0.0; 

Kvdot  =  0.0;  Kp       =  -1.10E-02*ZZ*LA2; 

Kr        =  -8.41E-04*ZZ*LA2;  Kv       =  3.055E-03*ZZ*L; 

%  YAW  HYDRODYNAMIC  COEFFICIENTS 

Npdot  =  0.0;  Nrdot  = -3.400E-03*ZZ*LA3; 

Nvdot  =  1.240E-03*ZZ*LA2;  Np       =  0.0; 

Nr       = -1.640E-02*ZZ*LA2;  Nv       =  -7.42E-03*ZZ*L; 

Zg  =  input(  enter  value  for  Zg') 

%  THIS  PERTURBATION  PROGRAM  SOLVES  FOR  THE  SOLUTION 
%  ABOUT  Xg=0 

Xg  =  0.0; 

index  =  1; 

for  Xg  = -1.0:0.05:0.2; 


a  =  (Ixx-Kpdot); 
d  =  (M*Zg); 


b  =  (Kp*U0); 
e  =  Kv*U0; 


g  =  (M*Zg*U0)+(Kr*U0);  h  =  Ypdot  4  <M*Zg); 

j  =  M  -  Yvdot;  k  =  Yv*U0; 

m=  (Yr*U0)  -  (M*U0);  o  =  Npdot; 

q  =  (Xb*BUOY)-(Xg*W);  r  =  (M*Xg)-Nvdot; 

w  =  (Nr*U0)-(M*Xg*U0);  x  =  Nv*U0; 


c=(Zg*WHZb*BUOY); 
f=  Krdot; 
i  =  Yp*U0; 
1  =  (M*Xg)  -  Yrdot; 
p  =  Np*U0; 
u  =  Izz  -  Nrdot; 


Al  =  (j*u)-(l*r);       Bl  =  -((u*k)+(j*w)-(m*r)-(l*i));      CI  =  (k*w)-(x*m); 

Ar  =  a;       Br  =  -b;       Cr  =  c; 

Kl  =  d*u*e  +  d*Nvdot*Kr*U0; 

K2  =  d*((Nr*Kv*U0A2KNv*Kr*U0A2)); 

alpha      =  -M*M*Zg*Kr*U0; 
beta        =  -d*((W*Yrdot)+(M*Kv*U0*U0)); 

gamma  =  (d*W*U0*(Yvdot-Yr))+(W*U0*(Kr*Yvdot-Kr*M-Kv*Yrdot)); 
delta       =  g*Yv*W*U0  -  m*Kv*W*U0; 
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%  P  Q  R  S  &T  SOLVE  THE  4TH  ORDER  POLYNOMIAL  FOR  aO 

P  =  Ar*Al; 

Q  =  Ar*Bl  +  Br*Al  +K1; 

R  =  Ar*Cl  +Cr*Al  +Br*Bl  +K2; 

S  =  Cr*Bl  +  Cl*Br; 

T  =  Cr*Cl; 

A  =  P; 

B  =  Q+alpha*Xg; 
C  =  R+beta*Xg; 
D  =  S+gamma*Xg; 
E  =  T+delta*Xg; 

z  =  [ABODE];  templ(index,l:4)  =  roots(z)'; 

xg(index,l)  =  Xg;  index  =  index  +  1;  end 

zz  =  [P  Q  R  S  T];  temp2(  1,1:4)  =  roots(zz)'; 

a0  =  temp2(l,4) 

%  THIS  PART  SOLVES  FOR  al  IN  THE  1ST  ORDER  PERTURBATION 

numl   =  -(alpha*(aOA3)+beta*(aOA2)+gamma*aO+delta); 
denl     =  (4*P*(aOA3)+3*Q*(aOA2)+2*R*(aO)+S); 
al  =  numl/denl 

%  THIS  PART  SOLVES  FOR  a2  IN  THE  2ND  ORDER  PERTURBATION 

num2a=-(alA2*(6*P*aO*aO+Q*aO+R)+al*(gamma+2*beta*aO+... 

3*alpha*a0*a0)); 
a2  =  num2a/denl 

index  =  1; 

for  X  = -0.8:0.04:0.4; 

%  PERTURBROOT1  IS   1ST  ORDER  PERTURBATION  APPROXIMATION 
perturbrootl(index)  =  aO  +  al*X; 

%  PERTURBROOT2  IS  2ND  ORDER  PERTURBATION  APPROXIMATION 
perturbroot2(index)  =  aO  +  al*X  +  a2*X*X; 

xlocation(index)  =  X; 
index  =  index+1; 
end 
return 
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